October 1996 
revised November 1997 
quant-ph/9712U52 



QUANTUM MECHANICS OF LATTICE GAS AUTOMATA 

11. BOUNDARY CONDITIONS 

: AND OTHER INHOMOGENEITIES 

a^ . 

o : 

CD ■ 

Q '. David A. Meyer 

m ; 

^ I Institute for Physical Sciences 

^ ! and 

CN . Project in Geometry and Physics 
^ ! Department of Mathematics 

CN ! University of Cahfornia/San Diego 
^ ; La Jolla, CA 92093-0112 

Q\ • dmeyer@chonji.ucsd.edu 

§ ■ ABSTRACT 

^3 : 

We continue our analysis of the physics of quantum lattice gas automata (QLGA). Pre- 
I vious work has been restricted to periodic or infinite lattices; simulation of more realistic 
^ I physical situations requires finite sizes and non-periodic boundary conditions. Further- 
ed ■ more, envisioning a QLGA as a nanoscale computer architecture motivates consideration 
of inhomogeneities in the 'substrate'; this translates into inhomogeneities in the local evo- 
lution rules. Concentrating on the one particle sector of the model, we determine the 
various boundary conditions and rule inhomogeneities which are consistent with unitary 
global evolution. We analyze the refiection of plane waves from boundaries, simulate wave 
packet refraction across inhomogeneities, and conclude by discussing the extension of these 
results to multiple particles. 
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1. Introduction 

Shor's discovery of a polynomial time quantum algorithm for factoring [1] stimulated a 
surge of interest in quantum computation (see the extensive bibliographies of [2]). Most 
work has concentrated on serial algorithms — sequences of unitary, few qubit* operations — 
the quantum version of serial Boolean logic [4]. Single quantum logic gates have been 
realized experimentally in ion traps [5] and quantum electrodynamics cavities [6], and 
short sequences of such unitary operations have recently been implemented with NMR 
[7]. All of these systems, as well as proposed solid state architectures such as arrays of 
quantum dots [8], exist physically in cZ > spatial dimensions and therefore naturally 
evolve in parallel. Imposing a single gate operation restricts the rest of the qubits to be 
invariant, i.e., they must evolve by the identity operator; at the opposite extreme all the 
qubits would evolve according to the same, local (few qubit) operation during a single 
timestep. A quantum computer evolving according to such a homogeneous, local, unitary 
rule would have the quantum version of the massively parallel architecture possessed, for 
example, by Margolus' CAM machines [9]. 

The simplest algorithms which would run on such an architecture are quantum cellular 
automata (QCA) [10] or quantum lattice gas automata (QLGA) [11]. Even in cZ = 1 
spatial dimensions QCA are capable of universal computation [12], and the existence of 
the universal reversible billiard ball computer [13] implies that QLGA are also, at least in 
d > 2 spatial dimensions. Just as classical LGA are most effectively deployed to simulate 
physical systems such as fluid flow [14] , however, QLGA most naturally simulate quantum 
physical systems [11,15,16]: with the simplest homogeneous evolution rule, one particle 
QLGA simulate the constant potential Dirac [11] or Schrodinger [17] equation, depending 
on the relative scaling of the lattice spacing and timestep. 

An earlier paper [15] initiated a project to analyze which physical processes QLGA 
can simulate effectively. In that paper and in this one we concentrate on the most general 
model for a single quantum particle with speed no more than 1 in lattice units, moving 
on a lattice in one dimension. The amplitudes for the particle to be (left, right) moving 
at a lattice point x & L are combined into a two component complex vector ?/)(t, x) := 
{^-i{t,x),'ip+i{t,x)) which evolves as 



■^(t +l,x) = W-i'lp{t, x-l) + wo'ip{t, x) + W+i'l/j{t, x + 1). (1.1) 



Here the weights Wi e M2(C) are 2x2 complex matrices constrained by the requirement 



* A qubit [3] is a quantum system whose state is a vector in a two dimensional Hilbert space, e.g., a 
spin-- particle fixed in space. 
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that the global evolution matrix 



Wo 



U := 



W-i 



W+i 
Wo 



(1.2) 



W-l 



w+l 



\ 



■■J 



be unitary. We showed in [11] that the most general parity invariant solution, up to unitary 
equivalence and an overall phase, is given by 



Describing the evolution by (1.1)-(1.3) assumes that the system is homogeneous in 
space and that the lattice L is isomorphic either to the integers Z or to a periodic quotient 
thereof, say Zyy. To simulate physical systems [18] more generally, the model should 
be extended to allow for finite size and non-periodic boundary conditions. Furthermore, 
envisioning a QLGA as a nanoscale quantum computer architecture [19,20,2] motivates 
consideration of inhomogcneities in the 'substrate', possibly as a step towards implementing 
logical gates [4] and away from simply simulating quantum physical systems. In [15] 
we showed how to introduce an inhomogeneous potential in the model; the purpose of 
this paper is to investigate more general inhomogcneities in the evolution rule, including 
boundary conditions. 

In Section 2 we consider the simplest possible modification of the evolution rule (1.1) 
for a boundary at a; = 0, say, setting W-i there to and allowing the weight wq to differ 
from the constant wq of the rest of the lattice. The resulting Type I boundary condition 
suggests the form for a corresponding Type I inhomogeneity where the global evolution 
matrix (1.2) is changed by replacing one of the wq blocks with a different matrix wq, and 
allowing the weights Wi to differ on either side of the antidiagonal through it. In Section 3 
we show that such an inhomogeneous rule is unitary provided 6 is the same for all the 
weights. 

There is a dual inhomogeneity across which p is constant but 9 may differ; we describe 
this Type II inhomogeneity in Section 4 and find the corresponding boundary condition. In 
Section 5 we observe that the Type I and II inhomogcneities can occur together, changing 
both p and 9. The corresponding Type III boundary condition has an extra degree of 
freedom, justifying distinct classification. 

In Section 6 we show how to find the eigenfunctions of U in the presence of these 
boundaries. In each case the result is a linear combination of left and right moving plane 




(1.3) 
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waves with the same frequency. On a finite lattice with two boundaries, the spectrum of 

U is discrete. In Section 7 we investigate the discrete spectra for pairs of each type of 
boundary condition, determining how they depend on the boundary parameters and what 
the consequences are for the eigenfunctions. 

Simulations of wave packets on lattices with boundary conditions and in the presence of 
inhomogeneities confirms that the physical consequences of these inhomogeneous evolution 
rules are as expected. We show some results in Section 8. 

We conclude in Section 9 with a summary and a discussion of the extension of this 
work to the multiple particle sector of the Hilbert space. 



2. Type I boundary conditions 



If our system is neither infinite nor periodic, we must model it on a bounded lattice, e.g., 
L = {xE'Z\0<x<N — 1}. Since there is no lattice point to the left of 0, it is clear that 
the evolution rule (1.1) must be adjusted there (as it must also be at the right boundary). 
Making the minimal change in the model, let us suppose that the global evolution matrix 
takes the form 

^ Wo W+i \ 

W-i Wq 



U := 



V 



(2.1) 



where the Wi are given by (1.3). Thus a left moving particle at x = 1 has the same 
amplitudes (given by Wj^i) to advect to x = and scatter to the left or right, and a right 
moving particle aX x — has the same amplitudes (given by W-i) to advect to x = 1 and 
scatter to the left or right, as each would were there no boundary. (The analogous form 
for the evolution rule at a right boundary is obtained by a parity transformation.) The 
only differences we allow for this Type I boundary condition are in the amplitudes for the 
evolution of a left moving particle at x = and for the scattering of a right moving particle 
at a; = which remains there during the advection step; these are given by Wq. 



The unitarity conditions UU"'' 



U^U impose the following constraints on wq: 



and 



= WqW^_]^ + lU+ilUo 

= w^^-JJJq + WqW_i. 



(2.2a) 
(2.26) 

(2.3a) 
(2.36) 



Let 



Wq 



yi y2 



(2.4) 
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Then, assuming cos p ^ 0, (2.26) implies 

y2 = -zcos6'sinp 
7/4 = sin 9 sin p, 

while (2.36) implies 

yi = iy3tajiO. (2.6) 
The normalization condition (2.2a) requires 

y3 = -ze*"cos^ (2.7) 

for some arbitrary phase angle e M. Combining (2.4)-(2.7), we find 

_ _ f e*"sin^ — zcos^sinp\ , . 

^° ~ l^-ze*^cos^ sin ^ sin p )' ^ ' 

which satisfies all the constraints (2.2) and (2.3). 

The Type I boundary condition defined by (2.1) and (2.8) gives the same amplitudes 
as (1.3) for the scattering of a right moving particle at x = which remains there; only 
the amplitudes for the scattering of a left moving particle at a: = differ from the no 
boundary situation. The latter depend on a single real parameter v characterizing the 
boundary. Notice also that these amplitudes do not vanish in the decoupled case p = 
(whence wq = 0). That is, Tuq 7^ is required to define unitary boundary conditions even 
when the particle has speed 1 everywhere else in the lattice. 



3. Type I inhomogeneities 

The boundary weight wq defined by (2.8) has the same form as the weight wq defined 
in (1.3), except that the factor of sinp in the first column is replaced by e*^. Thus we 
can interpret the evolution rule defined by (2.1) and (2.8) as describing a system where 
the coupling constant p satisfies cos p = at and to the left of x = 0. This would make 
W-i = = w+i, so there would be no advection to the left of a; = 0. This suggests that 
the Wq we found in Section 2 may be a special case of an inhomogeneity in the coupling 
constant p. So let us consider a Type I evolution rule inhomogeneity of the form: 



U 



v 



W_^ Wq 



(3.1) 



where the Wi =: Wi{p, 6) are defined by (1.3) and w[ := w^{p' , 9'). 
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Now the unitarity conditions impose constraints on the relation between the Wi and 
the as well as on the inhomogeneity matrix tug: 

= w'_^_-^^w'^_-^ (3.2a) 

= w'_-i^WQ + Wqw'^i (3.26) 

1 = w'_-^w'^-^ + WqWq + w_^_-^w^_^-^ (3.2c) 
= tUQtul]^ + w_^_iwl (S.2d) 

and 

= w'\w-^i (3.3a) 

= Wqw'_^_i + w'IiWq (3.36) 

/ = w^^w^i + w^Wq + w^_iW_i (3.3c) 

= w^.^^Wq + WqW_-^^. {3.3d) 



Constraint (3.2a) is automatically satisfied but, again assuming that cosp ^ 7^ cosp', 
(3.3a) requires sin(^ — 9') = 0, so we set 9' = 9. Using the form (2.4) for wq, we observe 
that the constraints {3.2d) and {3.3d) are the same as (2.26) and (2.36), so the j/j must 
satisfy (2.5) and (2.6). Constraint (3.26) requires that 

Vi = sin 9 sin p' 

, (3-4) 

j/3 = —i cos 9 sin p , 



which is consistent with (2.6), just as (2.5) is with (3.36). Combining (2.4), (2.5) and (3.4) 
we find 

^ ^ , I n \ ( sin ^ sin p' — zcos^sinp\ _^ 
^ Wq{p ,9 , p) := \ . a- , ■ n- h (3-5) 
^ ^ y— zcos^smp' sm^smp y ^ ' 

which also satisfies the remaining (normalization) constraints in (3.2) and (3.3). 

The arbitrary phase degree of freedom in the Type I boundary condition is not present 
in (3.5), but as anticipated, this Type I inhomogeneity describes a change in the coupling 
constant p, the mass 9 being held fixed across the inhomogeneity. The locus of the in- 
homogeneity is quite precise: a left moving particle from a; = obeys the 'primed' rules, 
while a right moving particle obeys the 'unprimed' ones. 



4. Type II inhomogeneities and boundary conditions 

The form (3.1) of the Type I inhomogeneity partitions the evolution matrix \J into two 
pieces across an antidiagonal through the wq block (inside the block the partition runs 
between the two columns). We might also consider an inhomogeneity which partitions U 
across an antidiagonal through a pair of W-\ and tu+i blocks. Such a Type 11 evolution 
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rule inhomogeneity has the form: 



U 



w 



w' 



w 



+1 



w 



-1 



Wn 



w 



+1 



w 



-1 



Wn 



w. 



+1 



v 



(4.1) 



where again Wi = Wi{p,$) and w'^ = Wi{p',$') are defined by (1.3) but with a priori 
different parameters. 

The unitarity conditions UU"'' = I = U^U impose even more constraints in this more 
comphcated situation: 






/ 



and 



/ = 






= 


w\_iWq + wlw_i 




= 


/f ^ 




= 






/ = 







(4.2a) 
(4.26) 
(4.2c) 
(4.2d) 
(4.2e) 

(4.3a) 
(4.36) 

(4.3c) 
(4.3rf) 
(4.3e) 



Suppose the inhomogeneity matrices have the most general forms: 

and tu+i := 



X3 X4 



(4.4) 



Then, assuming cosp 7^ 0, constraint (4.2c) requires Z2 = = Z4. Similarly, assuming 
cos p' 7^ 0, constraint (4.2(i) requires xi = = X3. Thus the inhomogeneity matrices have 
the same advection/scattering interpretation as in the homogeneous situation. 



Now constraints (4.3c) and (4.3d) imply that 

Z3 = izi tan 9' 



(4.5) 



and 
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X2 = ix^taxiO, (4-6) 

respectively. Imposing the normalization constraint (4.2a) we find that 

zi = e*^ cos p' cos 9'. (4.7) 

Then imposing the normalization constraint (4.3a) implies cos^ p = cos^ p', so we set 
p' = p. Combining (4.4), (4.5) and (4.7) gives 

w+^^e^^w'^,. (4.8) 
Similarly, imposing the normalization constraint (4.2e) we find that 

X4 = e*-*^ cos p cos 9. (4.9) 
Combining (4.4), (4.6) and (4.9) gives 

w-i = e''^w-i, (4.10) 

which also satisfies the last normalization constraint (4.3e). The two remaining constraints 
(4.26) and (4.36) require only that x = ~C (mod 27r), which can thence be set to by a 
unitary transformation. Thus (4.8) and (4.10) become 

W-i = W-i{p,9) and w+i = w+i{p,9'). (4-11) 

Just as the Type I inhomogeneity described by (3.1) and (3.5) specializes to a Type I 
boundary condition described by (2.1) and (2.8) when cosp' = so that there is no 
advection to the left of a; = 0, the Type II inhomogeneity described by (4.1) and (4.11) 
specializes to a boundary condition when cos 9' = 0. In this situation, when a left moving 
particle at x = 1 advects to .x = it scatters to the right, while a right moving particle at 
X — which remains at a; = also scatters to the right — a particle initially at x > or at 
a; = and right moving has no amplitude to be at x < or at x = and left moving at 
any subsequent timestep. 

This is a special case of the Type II boundary condition which we expect to be 
characterized by nontrivial phases, just as is the Type I boundary condition. The 'primed' 
parameters satisfy p' = p and cos^' = 0, so generalizing the Type II inhomogeneity by 
multiplicative phases suggests 



U := 



W-l Wo W+i 

W-i Wo 



(4.12) 



where 



_ , ie*^i sin^ . io\ 

«^-i=cosp( Q ) (4.13) 
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is a generalization of (4.10). Then the unitarity conditions impose constraints on the phase 
angles xi, X2-, v and C, via (4.2) and (4.3). Constraint (4.2c) is automatically satisfied, while 
(4.3(i) requires xi = X2 ='■ X (mod 27r). This means that the normalization constraints 
are necessarily satisfied so the only remaining constraints are (4.26) and (4.36). These are 
satisfied provided v = x C (mod 2tt). Finally, up to unitary equivalence we may set 
X = 0, so the most general Type II boundary condition is defined by: 



U = 







+1 



W-l 



Wq 
W-l 



W+1 
Wo 



\ 



■J 



(4.14) 



where = Wi{p,0). Just as does the Type I boundary condition, the Type II boundary 
condition has one phase degree of freedom. 



5. Type III boundary conditions 

The two types of inhomogeneities we have found reflect the p < — 6 duality evident in 
the dispersion relation (6.2) discussed in [15]: The Type I inhomogeneity has constant 9 
and discontinuity in p while the Type II inhomogeneity has constant p and discontinuity 
in 9. Suppose we wish to change both p and 9. This is clearly possible using a Type I 
inhomogeneity to change p followed by a Type II inhomogeneity to change 9, provided the 
discontinuities are sufficiently far apart that the constraints (3.2), (3.2), (4.2) and (4.3) 
do not overlap. In fact, the discontinuities can be adjacent: it is straightforward to verify 
that the evolution matrix 



U :-- 



Wr 



W 



-1 









Wo 


w+1 






Wo 






W-l 


Wo 



w 



+1 



(5.1) 



V ■••/ 

is unitary for W-i = W-i{p,9) and w+i = w+i{p,9') as in (4.11) and wo = wo{p' ,9' , p). 
The evolution matrix (5.1) describes a system in which the parameters p' and 9' change 
to p and 9 across the inhomogeneity. 

While (5.1) does not describe a new type of inhomogeneity as it is composed of a 
Type I and Type II pair, our experience with boundary conditions in the previous sections 
suggest that there may be an analogous Type III boundary condition which has extra 
phase degrees of freedom. Suppose 



U 



I Wo w+i 
W-l Wo w+1 
W-l Wo 

V 



(5.2) 
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where the Wi are given by (1.3), is given by (4.13), tu+i is generahzed from (4.11): 

_ / e^^i cos 9' 

W+i = cos P \ ■ iCo ■ Of n 

^ ^ \ie^^^ sm9 

and wq is the same as in (2.8) with 9 replaced by 9' and also with additional phase factors: 



Wq 



e*"i sin 9' — ■ie*"^ cos 9 sin p 
-ze*"3 cos 9' e*"* sin 9 sin p 



In this case, the unitarity conditions require 

/ = WqWq + Wj^-Jv^j^i (5.3a) 

= WqW^_^ + Wj^^Wq (5.36) 

= Wj^^w^_^ (5.3c) 

/ = w_-Ji]^_]^ + WqWq + Wj^-^w\^-^ (5.3ci) 

and 

/ = w\wq + w^_]W_^ (5.4a) 

= w'^j^^Wq + wyw_]^ (5.46) 

= w\{w_^ (5.4c) 

I — w\{Wj^-^ + WqWq + w'^_iW_i. (5.4(i) 

Constraints (5.3c) and (5.4c) are the same as (4.2c) and (4.36), respectively, so they have 
the same consequences as in the case of the Type II boundary condition: (5.3c) is satisfied 
automatically, while (5.4c) requires xi = X2 ='■ X (mod 2ti). Next, (5.36) implies V2 = 
X + ^i (mod 27r) and V4 = X + C2 (mod 27r). Constraint (5.46) requires v^ — vi = C2 — Ci —'■ ^ 
(mod 27r), whereupon the remaining (normalization) constraints in (5.3) and (5.4) are 
automatically satisfied. Combining these results and setting v := vi, C •= Ci gives 

,v f ism9\ _ ,r f cos 9' 



Wq 



w-i = e"^ COS p „ „ w+i = e"^ cos p . . „ 

VO COS 9 J ^ \ie^''sm9 

e^^ sin 9' cos 9' sin p 

.ieii^+S) cos 9' e^^^+^+^^ sin 9' sin p 



for the weights in (5.2). We may set x = = 5 by a unitary transformation, so the 
most general Type III boundary condition, up to unitary equivalence, is given by W-i = 
W-i{p, 9), tfJ+i = e*^t(;+i(p, 9'), and 

_ / e*" sin 9' — ie*^ cos 9' sin p 

Y —ie^^ cos 9' e*^ sin 9' sin p 

As expected, in addition to 9' there are two phase angle degrees of freedom: v and 
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6. Plane waves near a boundary 

The global evolution matrices (2.1), (4.14) and (5.2) describe unitary evolution of a single 
particle in the presence of a boundary of Type I, II, or III, respectively. Away from the 
boundary the local evolution is still given by (1.1) and (1.3), so the one particle plane 
waves 

we found in [15] still evolve, locally, by multiplication by e~^^^ at each time step, where cj 
satisfies the dispersion relation 

cos a; = cos A; cos ^ cos p + sin ^ sin p (6.2) 

and e e {±1}- In fact, any linear combination 

^Ij^''''\x) + A^P^-''^'\x) (6.3) 

evolves locally by phase multiplication as both k and —k satisfy (6.2) with the same 
frequency co. 

Consider the Type I boundary condition at a; = and suppose there is an eigenfunction 
iIj^'^\x) of the form (6.3), which should be interpreted as a linear combination of incident 
and reflected plane waves with relative amplitude A, just as in the situation of scattering 
off a potential step considered in [21]. Then 

wJoV'^'Ho) +m;+iV'^'"^(1) = e-^'^V^'^^O). (6.4) 
The linear combination (6.3) is well defined for a; < and 

w-iip^'^\-l) + WQiP^^\0) + w+iV'^'^n+l) = e-^^V'^'^Ho) (6-5) 
for any ^ e C, so subtracting (6.5) from (6.4) gives 

{wo - wq)iJj^^\0) = w-iijj^'^\-l). (6.6) 
Using (1.3), (2.8) and (6.3) in (6.6) we find 

^ _ (e*" - sin p) V'iV^ (0) - ^e"''' cos p V'+Y^ (0) 
(e^^ - sin p)V'i"i''''^ (0) - ie^'' cos p V'^Ti'''"^ (0) 

where 

^(.,e)^0) := ( . '''^f'Z^,^''~'''''T'"-^e.) (6-8) 

^ ^ y sm p sm ^ + e*'* cos pcosd — e J ^ ^ 

is the (unnormalized) eigenvector of D{k) in [15]. That is, with A given by (6.7), the linear 
combination (6.3) is an eigenfunction satisfying the Type I boundary condition. 
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The more comphcated Type II and III boundary conditions require modifications to 
the hnear combination of plane waves (6.3) near the boundary. Consider the Type II 
boundary condition and suppose 

V'^'^) (x) := V^'^'"^ {x) + ylV^"''''^ (x) for a; > 1 (6.9) 

and 

V'y(0):=0, (6.10) 

where the latter condition follows from the discussion preceding (4.12). At a; = 1 the same 
argument as in the Type I boundary case gives 

t(J_iV'^'^^ (0) = w-i {i;^''''^ (0) + ylV^"''''^ (0)) (6.11) 

which implies 

V'i^ (0) V'iV^(O) + M"i'''^(0)- (6-12) 
Applying (4.14) to the eigenfunction ip^^^x) at a; = gives 

e'^^w'o^^^^O) + e^^w;iV'^"^(l) = e'"'^ i^^'^^O) . (6.13) 

Using the expressions for w'^ with p' = p, cos^' = and (6.9), (6.10) and (6.12) in (6.13) 
we find 

^ _ (e^^ sinp - e-^'^)'i/'+i'^(0) + ie^('^+'=) cos p t(;^_!:f\o) 

{e< sinp - e-^'^)'i/'+"i^''Ho) + ie^^'^"^) cos pt(j^S['''^\o) ' 

Thus (6.9) with A given by (6.14), (6.10) and (6.12) define an eigenfunction satisfying the 
Type II boundary condition. 

To find the eigenfunctions for the Type III boundary condition, we still suppose they 
satisfy (6.9), but not (6.10). Since W-i is the same as for the Type II boundary, (6.11) 
still implies (6.12). Now applying (5.2) to the eigenfunction ip^^\x) at x = gives 

IfJoV'^'^nO) + w+ii;^'^\l) = e-*'^V^'^HO), 

which comprises a pair of linear equations for 'i/'^i''(0) ^ind A. These equations can be 
solved to give the eigenfunctions for the Type III boundary conditions, although we will 
not need the explicit solution here. 



7. Plane waves on finite lattices 

With only one boundary, e.^f., L = N as we were considering implicitly in the previous 
section, the wave number can take any value in the interval — tt < k < n and the fre- 
quency/energy spectrum is continuous with range 9 — p < < n — {9 + p) (assuming 
< p < ^ < 7r/2) determined by the dispersion relation (6.2). On finite lattices, however, 
the spectra are discrete and are determined by the two boundary conditions. Consider the 
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case of two Type I boundary conditions on a lattice of cardinality N. The weights in the 
boundary condition ai x = N — 1 are the parity transforms of those in (6.6): 

P{wo - wo)p-^ij^^\N - 1) = Pw-iP-^i;^^\N), 

where 

P:=(' 1 
\1 

This gives a second constraint on A: 

A = _^-2iKN-i) (e- - sin p) V^yr ^ (0) " i^^' cos p (0) 



^+i"\0) - ie''' cos pV'-i'^i 
(e»" - sinp)V'iri''''^(0) - ie-^'^ cosp vi"i^''^(0) 



which must be consistent with (6.7). To see how this determines the discrete spectrum, 
let = = p. Then (6.7) becomes 

A=-e — :t T : ^ (7.2) 

e-»'=cos6'-e-*<^'^ +sm6' ^ ^ 

and (7.1) becomes 

A = _e-2ik{N-i) cos g-e-^^'^ - sing 

e-^^ cos - e-*^'^ - sin 6* ■ ^ ' ' 

Setting the right hand sides of (7.2) and (7.3) to be equal and using the dispersion relation 
(6.2) to eliminate a;, we find, after some algebra, 

^-2i{N-2)k ^g-^ -is.in.k COS 0) ^sinO + i sin k cos 0. (7.4) 

Supposing k to be real, the right hand side of (7.4) is the complex conjugate of the paren- 
thesized expression on the left hand side, which implies that 

tan((iV-2)/c) +sinA;cot6' = 0. (7.5) 

The left hand side of (7.5) has poles at k = {n+ l/2)7r/ (AT — 2), n e Z, between each pair 
of which there must be a root of the equation. Thus (7.5) has N — 1 roots in the interval 
< A; < TT, giving A^— 1 distinct values for eigenfrequencies in the range 0—p < uj < tt — — p 
(assuming < p < < tt/2). But [/ is a 2A'' x 2N matrix so it must have 2N eigenvalues 
g-iw pigures 1 and 2 show the results of computing the eigenvalues of U numerically for 
N = 16: the eigenfrequencies are plotted as functions of the Type I boundary parameter 
V, set to the same value at each boundary. Notice that while most of the eigenfrequencies 
lie in the expected intervals, there are four which, over parts of the range of v, do not. 

To understand the origin of these unexpected eigenfrequencies, let us reconsider (7.4) 
and suppose that k has a nonzero imaginary part. Then for large N and the correct sign 
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Figure 1. The eigenfrequencies u! of U for a 
lattice of size N = 16 with p = 0, = tt/4 and 
two Type I boundary conditions with the same 
parameter v. 



Figure 2. The same situation as in Figure 1 but 
with parameters p = 7r/4, = 7r/3. In both cases 
there are two eigenvalues with w ~ when v — 
and two with u; ~ tt when v = tt. 



of k, the left hand side of (7.4) becomes arbitrarily small. So, if there were such a k which 
caused the right hand side of (7.4) to vanish, it would provide an additional root. Solving 

= sin 6' + -(e*^ - e"*'') cos6', 

we find 

e*'^ = -tan^isec^. (7.6) 

The negative root in (7.6) makes the norm |e~*'^| < 1 for < ^ < 7r/2; furthermore, 
it satisfies the dispersion relation (6.2) with iv — 0. Thus in the AT — > c>o limit, 1 is an 
eigenvalue of U with multiplicity two. For finite N these extra eigenfrequencies split, finely, 
and are only very close to 0. As wc sec in Figures 1 and 2, as f changes away from 0, 
the splitting increases and the eigenfrequencies move into the range associated with real 
wave numbers. An analogous discussion explains the pair of eigenfrequencies near tt at 
V = TT. The eigenfunctions having these eigenfrequencies corresponding to wave numbers 
with nonzero imaginary part are, of course, not plane waves; rather, each describes the 
state of a 'low' energy particle which is 'trapped' at the boundaries, with exponentially 
decreasing amplitude to be in the interior of the lattice. 

For the case of two Type II boundary conditions note that per the discussion following 
(4.12) the eigenfunctions of interest are those which have vanishing left (right) moving 
amplitude at the left (right) boundary. Thus when \L\ — N, there are 2A^ — 2 relevant 
eigenfunctions and eigenfrequencies. Figures 3 and 4 show the results of computing the 
eigenvalues of U numerically for A?" = 16: the eigenfrequencies are plotted as functions of 
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Figure 3. The eigenfrequencies lo of U for a 
lattice of size N = 16 with p = 0, 9 = n/A and 
two Type II boundary conditions with the same 
parameter 



CO 




Figure 4. The same situation as in Figure 3 but 
with parameters p = tt/4, = tt/S. In both 
cases there are two eigenvalues with w ~ and 
two with w ~ TT when = tt. 



the Type II boundary parameter (, set to the same value at each boundary. As in the 
Type I boundary situation, most of the eigenfrequencies he in the ranges corresponding 
to real wave numbers, although near C = there are four which do not, and which are 
explained by an analysis similar to that of the preceding paragraph. 

Finally, consider the case of two Type III boundary conditions, again with equal 
parameter values. In this case there is a non-phase parameter which can be adjusted, 
namely 9' . Figure 5 shows the eigenfrequencies of [/ as a function of 9' for the rule defined 
hy p — 0, 9 — IT /4:, with boundary parameters v — — To separate the eigenvalues 
we have computed them for a lattice of size only N = A. Figure 6 is similar, but the rule 
parameters are now p = w/A, 9 = tt/S. In this case N = 8 and the two eigenfrequencies 
near are only finely split over the whole parameter range. Notice that in each case 
there are actually six eigenvalues corresponding to imaginary wave numbers. Examination 
of the eigenfunctions shows that the two with eigenfrequencies near have amplitudes 
concentrated in the states |1, — 1) and \N — 2, +1), while the four with eigenfrequencies 
closer to ±7r have amplitudes concentrated at a; = and x = N — 1. 



8. Reflection and refraction of wave packets 

The physical meaning of the rule inhomogeneities we are considering is perhaps most clear 
in wave packet simulations. In [15] we defined binomial wavepackets with localized initial 
position and particularized initial wavenumber. In each of the simulations of this section 
the initial wavepacket is built from a plane wave (6.1) and (6.8) with ko = 7r/4, is centered 
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Figure 5. The eigenfrequencies uj of U for a 
lattice of size N = 4 with p = 0, = 7r/4 and 
two Type III boundary conditions with the same 
parameter O' and both boundary phase angles 0. 



Figure 6. The same situation as in Figure 5 
but with parameters p = 7r/4, 9 = tt/3, and 
N — 8. In both cases there are two eigenvalues 
with w ~ and four with a; ~ tt. 



at X = 16 and has width 32, on the lattice < x < 63. The peak frequency ujq and the 
group velocity depend on the rule parameters p and 9 through the dispersion relation (6.2). 

Let us first consider the reflection of such a wave packet from the possible boundaries. 
Figure 7 shows the evolution of the wave packet with parameters p — and 9 = n/A in the 
presence of Type I boundary conditions with v = 0. Reflection from Type II and Type III 
boundaries is extremely similar: in each case the signiflcant dispersion of the wave packet 
at the time of interaction with the wall results in a sequence of reflected (smaller) wave 
packets. 

As wc learned in [15], a 'massless' wave packet disperses more slowly than a massive 
one. In Figure 8 we show the results of a simulation of this case: p = n/A = 9 and 
the boundaries are both of Type II with C = 0. Reflection from Type I and Type III 
boundaries is again similar: in each case the wave packet reflects cleanly and suffers little 
more dispersion than if the boundary had not been there. 

Now let us consider the effect of rule inhomogenities on wave packet evolution. Figure 9 
shows the results of a simulation in which there is a Type I inhomogeneity at x = 31: the 
rule parameter 9 is constant at 7r/4 while p is to the left, and 7r/4 to the right, of 
the inhomogeneity. There is both reflection and transmission of the wave packet at the 
inhomogeneity: the reflected wave disperses rapidly which causes an interaction with the 
left boundary similar to that shown in Figure 7 while the transmitted wave packet has 
little dispersion and evolves much as the wave packet in Figure 8. 
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Figure 7. Evolution of the fcg = ''^/^ wave packet with width 32 for rule parameters p = 0, 
= 7r/4. The boundaries are both of Type I with v = 0. 



{figure available from author) 



Figure 8. Evolution of the same wave packet as in Figure 7 for rule parameters p = 7r/4 = 0. 
The boundaries are both of Type II with = 0. 
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Figure 9. Evolution of the same wave packet as in the previous figures with rule parameters 
= n/4 everywhere and p = to the left and p = 7r/4 to the right of a Type I inhomogeneity 
at a; = 31. Both boundaries are of Type I with v = 0. 



{figure available from author) 



Figure 10. Evolution of the same wave packet as in the previous figures with rule parameters 
p = 7r/4 everywhere and = 7r/4 to the left and = 7r/3 to the right of a Type II inhomogeneity 
at a; = 32. Both boundaries are of Type II with ^ = 0. 
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Next, let us consider the effect of a Type II inhomogeneity which changes the value 
of 9 from 7r/4 to the left of x = 32 to 7r/3 to the right. Figure 10 show the results of 
a simulation with p = n/A everywhere and Type II boundary conditions. To the left of 
the inhomogeneity the rule is 'massless'; this is evident in the negligible dispersion of the 
wave packet and its reflection off the inhomogeneity and then off the left boundary. With 
higher probability, however, the particle is transmitted through the inhomogeneity. The 
transmitted wave packet evolves according to a 'massive' rule and begins to disperse slowly 
so that reflection off the right boundary creates a small trailing wave packet. 

Recall from Section 5 that the Type I and Type II inhomogeneities can be adjacent 
with evolution matrix of the form (5.1). Figure 11 shows the results of a simulation of 
this situation: to the left of the inhomogeneity p = 0, 9 = 7t/4 and to the right p — 7r/4, 
9 — tt/3. The boundary conditions are of Type III with 9' — — v — (. There is 
the same concentration of probability at the inhomogeneity that occurs with the Type I 
inhomogeneity shown in Figure 9, together with less transmission and more dispersion to 
the right than with the Type II inhomogeneity shown in Figure 10. 

Finally, recall the p < — 9 duality displayed by the dispersion relation (6.2). Figure 12 
shows the results of a simulation in the presence of a Type III inhomogeneity constructed to 
convert the rule parameters p = tt/S, 9 = tt/4: on the left to the dual pair p = 7r/4, 9 = tt/3 
on the right. The reflected and transmitted wave packets have more even probabilities 
than in Figure 11 and evolve with the opposite group velocities. There is some asymmetry, 
most evident at the end of the simulation upon reflection from the boundaries; it is due 
to the asymmetry of the combined Type I/Type II inhomogeneity as well as of the initial 
condition. 

9. Discussion 

We have found dual inhomogeneities consistent with unitary global evolution of the general 
one particle rule (1.1)-(1.3): Type I implements a change in p while Type II implements 
a change in 9; adjacent Type I/Type II inhomogeneities implement changes in both p and 
9. Each of these three possibilities has a corresponding boundary condition character- 
ized by additional parameters. Despite this apparent variety of possibilities, the unitarity 
constraint is quite restrictive: besides the phase implementation of an inhoniogeneous po- 
tential [17,15] and some degenerate cases, these are the only possible local inhomogeneities 
up to unitary equivalence. 

A natural question to ask is how these rule inhomogeneities extend to the complete 
multiparticle rule set. Even the homogeneous rules for the general one dimensional QLGA 
whose one particle subspacc we have been investigating are quite complicated as there are 
effectively two, three and four particle interactions. In [11], however, we found the only 
particle preserving generalization of the rules (1.1)-(1.3) with p = (i.e., with particles 
of unit speed). In this case no more than two particles can advect simultaneously to a 
given lattice site, whereupon they scatter in opposite directions with amplitude e*"^, e R. 
In fact, the global evolution remains unitary if the constant phase angle considered in 
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Figure 11. Evolution of the same wave packet as in the previous figures with rule parameters 
p = 0, = 7r/4 to the left and p = 7r/4, 9 = 7r/3 to the right of a combined Type I/Type II 
inhomogeneity at a; = {Sl, 32}. Both boundaries are of Type III with = = v = C^. 



{figure available from author) 



Figure 12. Evolution of the same wave packet as in the previous figures with dual rule 
parameters p = n/3, = 7r/4to the left and p = 7r/4, 6 = 7r/3 to the right of a combined 
Type I/Type II inhomogeneity at a; = {31, 32}. Both boundaries are of Type III with 9 = = 
v = C 
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[22,11,17,21] becomes any function of the lattice sites ((>{x). That is, the two particle 
scattering amplitude can be any inhomogeneous function on the lattice, independently 
of the one particle 'scattering' amplitudes — and this independent inhomogeneity extends 
to the values of (f){x) at the boundaries. (Notice that even in the Type I and Type III 
boundary conditions where wq ^ 0, and hence a particle scattering off the boundary can 
have speed 0, no more than two particles can advect to the same lattice site.) The question 
of determining which boundary conditions are consistent with integrability of the model, 
via the Bcthe ansatz [23,24] as we began studying in [21] or by generalization of the Yang- 
Baxter equation [24] as has been used in closely related models [25], is of fundamental 
interest. 

For the purposes of quantum computation with QLGA, we conclude by noting that we 
have explicitly formulated the possible local inhomogeneities in the one dimensional unit 
speed model. Extending these multiparticle results to multiple speeds and higher dimen- 
sions seems likely to be algebraically more complicated but conceptually similar — single 
particle single speed models with inhomogeneities have been constructed in two dimen- 
sions [26,27,17]. More interesting is the question of how to exploit such inhomogeneities to 
effect specific quantum computational tasks more efficiently than by simply implementing 
a quantum version of reversible billiard computing [13,20] using a homogeneous rule. The 
most natural use of QLGA may be to simulate other quantum physical systems; designing 
an inhomogeneous QLGA to be an efficient universal quantum computer may consequently 
be difficult. A reasonable intermediate goal would be to solve specific problems particu- 
larly well suited to this architecture. Although neither implements a quantum algorithm, 
Squier and Steiglitz' particle model for parallel arithmetic [28] and Benjamin and John- 
son's recent proposal for an inhomogeneous nanoscale cellular automaton adder [29] may 
provide useful points of departure. 
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